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Abstract 

We study equilibrium time correlations for the discrete nonlinear Schrodinger 
equation on a one-dimensional lattice and unravel three dynamical regimes. There is 
a high temperature regime with density and energy as the only two conserved fields. 
Their correlations have zero velocity and spread diffusively. In the low temperature 
regime umklapp processes are rare with the consequence that phase differences appear 
as an additional (almost) conserved field. In an approximation where all umklapp is 
suppressed, while the equilibrium state remains untouched, one arrives at an anhar- 
monic chain. Using the method of nonlinear fluctuating hydrodynamics we establish 
that the DNLS equilibrium time correlations have the same signature as a generic 
anharmonic chain, in particular KPZ broadening for the sound peaks and Levy 5/3 
broadening for the heat peak. In the, so far not sharply defined, ultra-low temper¬ 
ature regime the integrability of the dynamics becomes visible. As an illustration 
we simulate the completely integrable Ablowitz-Ladik model and confirm ballistic 
broadening of the time correlations. 
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1 Introduction 


The one-dimensional discrete nonlinear Schrodinger equation (DNLS) is widely used as 
theoretical description of nonlinear optical wave guides [1], semiclassical approximation to 
the dynamics of Bose-Einstein condensates [2] , electronic transport in biomolecules [3] and 
more, as discussed in the surveys [4, 5]. The DNLS has a surprisingly rich dynamical behav¬ 
ior. Most studies are concerned with finite energy solutions, which physically corresponds 
to zero temperature. In this article we focus on the dynamics at non-zero temperatures, a 
topic which lies in the domain of non-equilibrium statistical mechanics. 

DNLS governs the dynamics of a complex-valued wave field denoted by ifjj with j 6 Z, 
the one-dimensional lattice, and is defined through the hamiltionian 

^ = ZX^i^+ i -^i 2+ ^i^i 4 ) • ( L1 ) 

jez 

We will consider the defocusing case with g > 0. Kulkarni and Lamacraft [6] recently 
simulated the DNLS dynamics at density p = 1, m = 1, g = 1, and temperature j3~ l = 
0.005. They discovered that the density-density time correlations in thermal equilibrium 
have symmetrically located sound peaks, which travel ballistically and broaden as t 2 / 3 , 
consistent with the predictions based on the KPZ (alias stochastic Burgers) equation. As 
further amplified in [7], one writes down a pair of coupled noisy Burgers equations for the 
left and right moving modes. Since the peaks separate linearly in time, it is argued that 
they decouple and each peak is governed by a scalar noisy Burgers equation, for which 
the exact scaling exponent and shape function is known [8, 9, 10, 11]. In contrast, the 
high temperature DNLS has diffusive transport, as confirmed through molecular dynamics 
(MD) simulations of a steady state with stochastic boundary conditions forcing a flux of 
density and energy [12]. Measured are the Onsager coefficients relating the linear response 
in the flux to the forcing. As we will argue, there is no sharp transition but the borderline 
between low and high temperatures is roughly characterized by 

Pgp 2 ~ 4, 1- (1-2) 

Increasing the temperature even further one reaches the infinite temperature line, in de¬ 
pendence on the density, which borders the region of negative temperature states, realized 
through the appropriate microcanonical ensemble. There one observes interesting coarsen¬ 
ing processes which are carried by the motion and coalescence of breathers [13]. 

In our contribution we will develop a global picture for positive temperature states, more 
tightly linked to DNLS than the previous discussions [6, 7]. The central dynamical concept 
are umklapp processes, at which the phase difference between neighboring sites crosses the 
values ±7r. At high temperatures umklapp happens frequently. Density and energy are 
the only conserved fields and their correlations spread diffusively with no systematic drift. 
However, in the equilibrium ensemble at low temperatures the phase differences are small, 
of the order l/\//L Umklapp is an activated dynamical process and occurs only with a 
frequency of the order e -/3AE , with A E a suitable energy barrier. As the temperature is 
lowered, the field of phase differences becomes almost conserved. The conservation law for 
phase difference degrades as e~ t//r with a life time r proportional to e^ AE . For all practical 
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purposes, in particular for MD simulations, the phase differences are locally conserved, once 
f3 is in the domain with a border line specified by (1.2). Thus one arrives at a system of 
three local conservation laws. 

The described dynamical mechanism is familiar from second sound in liquid helium. 
There is no first sound for DNLS. But sound propagation is enabled through the appearance 
of an additional conserved field. In one spatial dimension such a transition is not sharp. 

In fact there is a third regime, baptized ultra-low, whose precise border has still to 
be explored. In the low temperature regime there are three sharp peaks and for them we 
expect universal scaling laws with the dependence on the microscopic interactions lumped 
together in a few non-universal coefficients. In the ultra-low regime the integrability of the 
dynamics becomes visible. The density-density time correlations have now in addition, or 
instead, broader structures which self-similarly expand on a scale linear in time. 

For the low temperature regime we elucidate the universal features, using the method 
of nonlinear fluctuating hydrodynamics. For this purpose our main tool is an effective 
low temperature hamiltonian, build such that the equilibrium ensemble remains unchanged 
while umklapp is completely suppressed. Hence phase difference, density, and energy are 
strictly conserved. In the well understood case of an anharmonic chain [14], there are also 
three conserved fields, stretch, momentum, and energy. The equilibrium state is of product 
form in the lattice index j. The DNLS hamiltonian contains however a nearest neighbor 
interaction, which complicates considerably the computation of the coupling constants for 
the normal mode representation used in nonlinear fluctuating hydrodynamics. Still we 
manage to ensure that the self-coupling of the heat mode vanishes, while the corresponding 
coefficients for the sound modes are different from zero. With this input the spreading and 
shape of the three peaks are predicted to be identical to the ones of a generic anharmonic 
chain, up to model-dependent non-universal coefficients. 

In essence, our predictions are based on conservation laws and a sufficiently chaotic 
nonintegrable dynamics. With this view, one would expect that our results are also appli¬ 
cable to one-dimensional quantum systems. In fact, quite some time ago Andreev [15, 16] 
studied low temperature one-dimensional Fermi fluids and argued already that the sound 
peaks broaden as f 2//3 , which turns out to be the same scaling exponent as for DNLS. We 
refer to [17] for more recent studies. One-dimensional quantum fluids in the continuum 
are not readily accessible to quantum simulations. On the other hand, by quantizing the 
DNLS hamiltonian one arrives at the Bose-Hubbard model on a one-dimensional lattice. 
The propagating sound peaks are expected to be visible at low temperatures and density 
of order one, which are parameters more favorable for quantum simulations. Of course, 
more optimistically one would hope to be able to realize such propagating sound peaks in 
an cold atom experiment. 

In view of the ultra-low temperature regime, we simulate the Ablowitz-Ladik model, 
which is an innocent modification of the DNLS dynamics, but known to be completely 
integrable. For this model, in a certain sense, every conserved mode generates a peak. 
Thus heuristically the correlations should be extended and widen ballistically. For the 
continuum NLS the same time correlation structure is expected to appear. The underlying 
lattice is essential for our findings. 

To provide a brief summary: First the high temperature regime is studied and simulated 
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at [3 — 1, with all other parameters as before. We then transform to action-angle variables 
and determine the low temperature Hamiltonian. For it we compute the nonlinear coupling 
coefficients in the normal mode representation and determine the Landau-Placzek ratios. 
Thereby the predictions for the shape functions and the non-universal coefficients are made 
available, which are then compared with MD simulations of DNLS at /3 = 15. For the ultra- 
low temperature regime we set f3 = 200. In a final chapter we discuss the Ablowitz-Ladik 
model. 


2 One-dimensional DNLS, basic properties 

Our starting point is the complex-valued field ipj governed by the DNLS hamiltonian 

N-l 

H = + U\^j\ 4 ) ( 2 - 1 ) 

l=o 

with mass m > 0 and coupling constant g > 0, i.e., a defocusing nonlinearity. We impose 
periodic boundary conditions, ?/;,v = V’o• Our interest will be in random initial data dis¬ 
tributed according to a thermal Gibbs state, for which we will take the limit N —> oo at a 
suitable stage. The dynamics is defined through 

i ah’.f = V f ’ < 2 ' 2 > 

where * denotes complex conjugate. Then 

i = “2+ 9 Hj \' Vh ( 2 -3) 

with the lattice Laplacian A = —d T d and dipj = ijjj. The kinetic energy is chosen such 

that in the limit of zero lattice spacing one arrives at the continuum nonlinear Schodinger 
equation. One could also expand the square, resulting in a hopping term plus a contribution 
proportional to the particle number, 


N-l 

N = (2.4) 

3=0 

The sign of the hopping term does not play a role, since it can be flipped through the gauge 
transformation tpj i—> e 17rj '0j. 

The DNLS has two obvious locally conserved fields, density and energy, 

Pi = iVhf , e j = ~ Vhf + l 9 I ipj I 4 • (2-5) 

According to the discussion in [18], the DNLS is nonintegrable and one expects density and 
energy to be the only locally conserved fields. They satisfy the conservation laws 

d iPj^f) + J P ,j+i(t) — = 0 , (2 6) 

^~ t e j{t) + J e j + i(t) — = 0 , 
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with density current 


Jm = ytwj-i i -'s; 3^-1) 


(2.7) 


and energy current 

Jej = ^j-1 - ■ (2.8) 

Canonically conjugate variables are introduced by splitting the wave function into its 
real and imaginary part as 

+ Wj) ■ ( 2 - 9 ) 

In these variables, the hamiltonian reads 


V-l 

H = J2 (i^((%) 2 + ( d Pj ) 2 ) +1 g{q] + p 2 j) 2 ) ■ ( 2 -io) 

j= 0 

The dynamics governed by Eq. (2.3) is then identical to the hamiltonian system 

fai = e n n, t-, Pl = ~a„H. ( 2 . 11 ) 

Note that H is symmetric under the interchange qj <H- p r 

It will be convenient to make a canonical (symplectic) change of variables to polar 
coordinates as 

<Pj = arctan (pj/qj), pj = | (p) + qf) , (2.12) 

which is equivalent to the representation 

V’i = • (2.13) 

In the new variables the phase space becomes (pj,Pj) € M+ x S 1 , with S 1 the unit circle. 
The corresponding hamiltonian is given by 

N-l 

H = Y1 ( 2 k(VPi+i Pj 2 (! - cos(^- + i - pj)) + (v ^+1 - Vp]) 2 ) + \9p ]) (2.14) 

j =0 

N-l 

= cos(p j+ 1 - Pj) + ± Pj + \ gpf) . (2.15) 

j= 0 

The equations of motion read then 

= & Pi = d Vi H. (2.16) 

From the continuity of yq(f) when moving through the origin, one concludes that at p 3 {t) = 
0 the phase jumps from Pj(t) to Pj(t) + tt. The pj's are angles and therefore position-like 
variables, while the pj s are actions and hence momentum-like variables. The hamiltonian 
depends only on phase differences which implies the invariance under the global shift pj 1 —> 

Pj + 0. 
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3 Equilibrium time correlations at high temperatures 

We consider the DNLS in thermal equilibrium, as described by the canonical ensemble 

N/2—1 

Z N (n,/3)~ 1 e~^ {H ~ flN) JJ (hlJjd'i//* (3.1) 

j=-N/2 

with inverse temperature f3 and chemical potential /r, /3 > 0, /i G M. For the theory we take 
the limit N —> oo. Numerically N will be set to 4096. Some of the derivations are done 
first at finite volume. Notationally we do not distinguish between finite N and N = oo, 
assuming that it will be understood from the context. Our interest is the propagation of 
small perturbations at large space-time scales, which are encoded by time correlations of 
the conserved fields, in our case density and energy, 

Sppti, t) = ( pj(t ); p 0 (0)), SfJJj, t ) = (pj(t); e 0 (0)), S ee (j, t ) = (e^t); e o (0)). (3.2) 

Here (•) denotes the average and (•; •) the second cumulant with respect to the canonical 
state (3.1). S(j,t) is defined on the entire lattice Z. By invariance under reversal of time 
and space, S pe (j,t ) = S ep (j,t). 

The canonical ensemble is even under complex conjugation of 'i/jj , while the density and 
energy currents are of the form i (z — z*). Hence 

(J P j) = 0, {J e j)= 0. (3.3) 

If the dynamics is sufficiently chaotic, this property signals diffusive transport. More pre¬ 
cisely we denote by u p the random deviation of the density from the uniform equilibrium 
density and by u e the one for the energy density. Then linear fluctuating hydrodynamics 
asserts that the random deviations are governed by the linear Langevin equation 

d t u a + d x { - d x (Du) a - (B£) a ) = 0, (3.4) 

a — p,e. Here D is the 2x2 diffusion matrix, B the noise strength matrix, and £ a (x,t) 
are two independent normalized space-time white noises. We define the static covariance 
matrix through 

OO 

C ao > = 5 ««'(.?,0), (3-5) 

j=-oo 

from which it follows that C = C T , C > 0. The assumed fluctuation dissipation relation 
reads then 

2DC = BB t , (3.6) 

which implies that the Onsager matrix DC is a symmetric, strictly positive matrix. In 
particular D has strictly positive eigenvalues. The covariance of the stationary process for 
(3.4) reads 

K(MK,(0,0)) = s aa .(x,t) = [ d ke- i2 * kx (e-^ 2Dt C ) aa ,. (3.7) 

J R 
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Figure 1: Equilibrium time-correlations S aa '(j,t ) of the density and energy at high tem¬ 
perature, /3 — 1, for the discrete NLS with initial states drawn from the canonical ensemble 
(3.1). The black dashed curves are the entries on the right of Eq. (3.8). 


Working out the Fourier transform, one arrives at the prediction 

< 3 - 8 ) 

The square root of D is unambiguous, since D has strictly positive eigenvalues. 

In Fig. 1 we show the results of a MD simulation for N = 4096, m — 1, g — 1, /3 — 1, 
and ( pj) = 1. The black dashed curves show the right side of (3.8), with C determined by 
the sum rule (3.5) and D fitted numerically, under the constraint of DC being symmetric. 
The numerical entries of C are 


/0.580 0.907\ 
\0.907 1.848J ’ 


(3.9) 


and we record D in the following table: 

t 256 512 1024 1536 

— / 3.145 —0.353\ / 3.137 -0.364\ /3.103 -0.356\ 73.079 -0.350\ 

^2.509 0.822 ) ^2.415 0.853 ) ^2.350 0.876 ) ^2.298 0.897 ) 

The values for D drift only little in time, which supports the diffusive scaling exponent 
consistent with linear fluctuating hydrodynamics. In Fig. 1 only the largest time, t = 1536, 
is shown. 


4 Dynamics at low temperatures 

In the MD simulations [6] the parameters of the hamiltonian are the same as used here, but 
/3 = 200 at system size N = 16384. For the density-density correlations a behavior very 
different from the one reported in Sect. 3 is observed, in exhibiting ballistic propagation 
and non-diffusive spreading. A main goal of our contribution is to understand and predict 
such low temperature properties on the basis of nonlinear fluctuating hydrodynamics. The 
argument is involved and we choose to first illustrate our method for the case of coupled 
rotators (CR), for which the dependence on the analogue of pj is much simpler. 
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Coupled rotators. The CR hamiltonian reads 


N -1 

Hck = Y (b? - cos(^- +1 - (fj)) . (4.1) 

j =o 

Here pj G S' 1 is the angle and Pj G R the angular momentum of rotator j. We impose peri¬ 
odic boundary conditions as p^ = Po- When compared to (2.14), the prefactor p^Jpj+i Pj 
has been set equal to 1 and the constraint pj > 0 is ignored. Angular momentum and en¬ 
ergy are the only locally conserved fields for Hqr- Under the canonical equilibrium measure 
their average currents vanish as N —» oo, signaling diffusive transport at high temperatures, 
as has been confirmed through MD simulations [19, 20]. 

At zero temperature, there is the one-parameter family of ground states with pj = p, 
Pj = 0. When heating up, under the canonical equilibrium measure, the phase pj jumps 
to Pj+i with a jump size 0( l/\/]d). For a more quantitative version, let us introduce the 
phase difference fj = 0(pj +1 — ipj), where 0 is 27r-periodic and Q(x) = x for |x| < 7r. 
Since 0 has a jump discontinuity, fj is not conserved. In a more pictorial language, the 
event that \pj + \(t) — Pj(t )| = 7r is called an umklapp for phase difference fj or an umklapp 
process to emphasize its dynamical character. At low temperatures a jump of size 7r has a 
small probability of order e~P AV with AV = 2 the height of the potential barrier. Hence 
fj is locally conserved up to umklapp processes occurring with a very small frequency only, 
see [19] for a numerical validation. 

In the low temperature regime it is tempting to use an approximation, where the po¬ 
tential V(x) = — cosx is Taylor expanded at the minimum x = 0. But such procedure 
would underestimate the regime of low temperatures, as can be seen from the example 
of a potential, still with AV = 2, but several shallow minima. The proper small param¬ 
eter is /3 _1 such that (3AV > 2, where 2 is chosen to have a safety margin. To arrive 
at an optimal low temperature hamiltonian, we first parametrize the angles po,... ,Pn-i 
through Tj = <Pj+i — Pj with rj G [—7r,7r]. To distinguish, we denote the angles in this 
particular parametrization by 4>j . The dynamics governed by Hqr corresponds to periodic 
boundary conditions at rj = For a low temperature description we impose instead 
specular reflection, i.e., if Tj = ±7r, then pj, p 3+ \ are scattered to p) = Pj+i, p'j +1 = Pj . By 
fiat all umklapp processes are now suppressed, while between two umklapp events the CR 
dynamics and the low temperature dynamics are identical. The corresponding hamiltonian 
reads 

2V-1 

ffcwt = Y. (Id + U*+1 - <h)) ( 4 . 2 ) 

3=0 

with 

V(x) = — cosx for |x|<7r, V(x) = oo for |x|>7r. (4.3) 

As before, periodic boundary conditions (pN = are understood. The pair ( (pj,Pj) are 
canonically conjugate variables. Note that as weights exp[— /3H CR ] = exp[—/3i?cR,it]- Thus 
all equilibrium properties of the coupled rotators remain untouched. 

The hamiltonian Hcr,h is a variant of the hard collision model with square well potential, 
see [21], thus a conventional anharmonic chain. The dynamics governed by idcR,it has 
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three conserved fields, the stretch r 3 = (pj+± — (pj , the momentum pj and the energy e 3 = 
\p 2 j + V{'f-j ). To study its equilibrium time correlations the methods and results from [14] 
apply directly. Because of (p 0 = <Pn, one has r j = 0 - The model is in the dynamical 

phase characterized by an even potential at zero pressure. 

We claim that, for j3 AU > 2, the CR equilibrium time correlations are well approx¬ 
imated by those of -£/cR,it, provided the time of comparison is not too long. The latter 
correlations can be obtained within the framework of nonlinear fluctuating hydrodynamics. 
Thereby one arrives at fairly explicit dynamical predictions for the low temperature regime 
of the CR model. 

Low temperature approximation for DNLS. We return to the DNLS hamiltonian 
and restrict ourselves to the case p > 0. As can be inferred from the hamiltonian (2.14), in 
the limit (5 —> oo the canonical measure converges to the one-parameter family of ground 
states with pj = p = p/g, (fj = p and Cp uniformly distributed on S 1 . As for CR, we 
introduce fj = @(<pj + i — tpj). In the context of Bose-Einstein condensates v 3 = —fj is 
called the superfluid velocity. At low temperatures, pj — p is confined by an essentially 
harmonic potential, hence \pj — p\ — (D{l/\/]3) and the phase p 3 jumps to <Pj+i with a 
jump size 0( l/\//3), also. 

The low temperature hamiltonian is constructed in such a way that the equilibrium 
ensemble remains unchanged while all umklapp processes are suppressed. To achieve our 
goal we follow verbatim the CR blueprint. Now (< pj,Pj) are a pair of canonically conjugate 
variables, only p 3 > 0 instead of p 3 e M. The phases are parametrized such that <Pj+i — Pj 
lies in the interval [—7r, 7r]. Umklapp is a point at the boundary of this interval. Thus the 
proper low temperature hamiltonian reads 


N—l 




{VPj +1 Pj U (^+1 - (Pi) + V (pj)) > 

j-0 

(4.4) 

where 





U(x) 

= —Acos(x) for \x \< tt , U(x) = oo for \x \> n , 

(4.5) 

and 

V(x, 

) = ^ xJ r\gx 2 for x>0, U(x) = oo for x<0. 

(4.6) 


For some computations it will be convenient to replace the hard collision potentials U, V 
by a smooth variant, for which the infinite step is replaced by a rapidly diverging smooth 
potential. The dynamics is governed by 


= = (4.7) 

including the specular reflection of p 3 at p 3 = 0 and of rj at r 3 = ±7r. As for CR, between 
two umklapp events the true and the low temperature dynamics are identical. Also as 
weights exp[— j3H] = exp[—/3-ffit], which is a salient feature of our approximation. 

There are two potential barriers, A U and AU. The minimum of V(x) — —x — px is at 
p — p/g , hence AU = \gp 2 . The minimum of U is at <f>j+\ — (pj — 0 and, setting p 3 = p, 
one arrives at A U = Ap. Thus the low temperature regime is characterized by 

|/3p 2 >2, lBp>2. (4.8) 
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In this parameter regime we expect the equilibrium time correlations based on H\ t to well 
approximate the time correlations of the exact DNLS. 

As for a conventional anharmonic chain we introduce the stretch rj = 4>j+\ — (j)j- Then 
the conserved fields are pj, fj , and the energy 

e i = v^i+i Pi U ( r i) + v ( Pi ) • ( 4 - 9 ) 

This local energy differs from the one introduced in (2.5) by the term ^(pj +1 — Pj)- hi 
the expressions below such a difference term drops out and in the final result we could use 
either one. The local conservation laws and their currents read, for the density 


with local density current 

dt.Pi + ^pj+i ^pP ^ 

(4.10) 

for the stretch 

<^pd = VPi- 1 Pi U ( r i-i) > 

(4.11) 

with local stretch current 

d t r i 3r,j- 1-1 — <Jr,j = 0 

(4.12) 

J r j 2 \j Pi+^l Pi ^( r i) T 2 yj Pi- 1/ Pi U( r i- i) + V (Pj ) , 

and for the energy 

(4.13) 

with local energy current 

dt e i + i7e,j+l <Je,j 0 

(4.14) 

— 2 VPi- 1 Pi+ 1 {U( r i- 

-i)V(r s ) + U\rj-i)U(rj)) + VftUft Ujr^Vjn). 

(4.15) 


To shorten notation, we set gj = (pj,rj,ej) and J 3 = (J'pj, J r ,j, J e ,j) ■ 

To stress again, the low temperature hamiltonian H\ t serves only as an input to nonlinear 
fluctuating dynamics. Our MD simulations will use DNLS throughout. 


5 The coupling coefficients for fluctuating hydrody¬ 
namics 

To arrive at predictions for the dynamics based on H\ t we will adopt the method from [14]. 
In contrast to a standard anharmonic chain, the hamiltonian, written in terms of stretches, 
is no longer a sum of one-particle terms. At first sight this seems to complicate matters 
considerably. But, there is a surprising identity for the three average currents, which allows 
us to still obtain the nonlinear coupling coefficients in a concise form. 
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The case of general p, is. The chemical potential /i is the dual variable for pj and we 
introduce is as the dual variable for Vj. Since </> 0 = 4>n, one has r j = 0- Hence, by the 

equivalence of ensembles, the equilibrium state is actually at v = 0. But for the quadratic 
expansion of the Euler currents it is more convenient to first work with general is, setting 
is = 0 at the end. The canonical ensemble of H\ t for a finite system with N lattice sites 
and periodic boundary conditions is given by 

N-l 

Z N (ji, u, (3)- 1 Po-^U 1 rj) dp . dr . (5.i) 

3 =0 

with the normalizing partition function 

r N-l 

Z N {p, u,P)= / e -/3(H- M EfTo 1 Pj-v EfTo 1 rj) TT dp . dr . ( 5 . 2 ) 

J (R+Xl-TT,*])" j= q 

On purpose, the index “It” has been omitted, since H and H\ t define the same measure. 
The canonical free energy is defined as 

F(/i, is, 0) = -/T 1 lim A log Z N (p, is , p ). (5.3) 

V—>-oo 

The averages of pj, rj, ej are 


P = (Pj) = ~ d n F (P, 0/3), r = (rj) = -d v F(p, is, p), 
e = i im n( H ) = d p(P F (h, O /3)) + p p + is r, 

iv—>-oo 


independent of j by translation invariance. By convexity of F , these relations define the 
inverse mapping (p, r, e) i—> (p(p, r, e), is(p, r, e), /3(p, r, e)). As discussed in Appendix B, the 
average currents turn out to be 

((>-fp,ji ffr,j , ffe,j)) (O P'-i p v) j. ('J-’j) 

The product form of the energy current will be crucial in what follows. 

Taking derivatives with respect to p, is, j3 generates sums over j. We therefore introduce 
the shorthand 

N-l 

(ifo’, h)) = ( 5 - 6 ) 

3=0 

with (•; •) denoting the second cumulant. Here / 0 refers to some local function and fj 
is the same function shifted by j. At the end of the computation one wants to take 
N —» oo. Infinite volume correlation functions, such as ho) = (fj ho) — (fj)(ho), decay 
exponentially fast to zero. Hence all our formulas hold also for infinite volume. Note that 
in this case 

OO 

<(/o;fto)>= £ </>;*»>- (5-7) 

j=-oo 
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With such notation the static susceptibility matrix is defined through 

/<(Po; Po)> <(po; ro)) ((po; e 0 ))\ 

Caa' = ((g a , 0 \9a', 0 » , C= \ ((^h Po)) (4 05 r o)) (4o; &o)) , (5.8) 

\ ((eo; Po )) (( e o; r o)> ((eo; e 0 )) / 

a, ol = 1, 2, 3, where the vector g 3 of the field variables is defined below Eq. (4.15). 

The following relations hold, 


In matrix notation, 


'd,, 


,4 


= ~((fo] eo — p po — v r 0 )). 

(5.9) 

((/o; Po))\ 


<</o; t- 0 >> ■ 

(5.10) 

((/o;e 0 ))/ 



Hence, using (5.5), one deduces that 


Setting 


0 1 v 

((57«,0) 9a', o)) W I 1 0 P 

^ ' V yU 2 /jLU, 


N -1 iV-1 

<(/o; 4; 4)) = E (/o; 4; 4), «(/o; 4; 4))) = ^ to; 4; 4) 

1=0 *,l=o 


(5.11) 


(5.12) 


for the third cumulant, one obtains 


<4(/o; 4) = /4(/o; 4; po)), 4(/o; M = /5<(/o; 4; r 0 )), 

4) = -((/o; 4; eo - PPo - vr Q )). 

Employing (5.11), one arrives at 


hi P al — (((i7l,o; 9a, 0 ] <7 7 ,o))) 


hi r ai — (((572,o; Po.o; p 7 ,o))) 


%a 7 — (((573,0; Pa, 0 ; P 7 ,o))) 



(5.13) 


(5.14) 

(5.15) 

(5.16) 


To write down the equations of nonlinear fluctuating hydrodynamics, one needs the 
average currents (5.5) expanded to second order at the uniform background values (p, r, e). 
To streamline our formulas, we denote the deviation from the background as u — (ui,U 2 , M 3 ) 
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= (p — p, r — r, e — e). The background values will be suppressed in our notation. To linear 
order, 

/ d p v d r is d e is \ 

^oa 1 ^u a i]a ( dpp, d r p d e p I • (5.17) 

\d p (pis) d r {pu) d e (pis)J 


We use the identity 

(dpp d p x d fl e\ / d p p d p v d p 0\ 

J d v p d v x d u e J d r p d r v d r p = t (5.18) 

\d,gp dp x d 6 e) \d e p d e is d e /3J 

to express derivatives of p, is, f3 by derivatives of p, r, e. Applying to the left matrix the 

relation (5.10) leads to 


(dpp <9 M r d lt e\ fP 0 0 \ 

d v p d v x d u e = 0 P 0 C. (5.19) 

\dpp dpx dpej \p is -1/ 


By matrix inversion, one thus obtains the right matrix in Eq. (5.18) as 


fdpp dpis d p p\ 1 [l 0 0 \ 

J d r p d r is d r P = — C I 0 1 0 J . (5.20) 

\d e p d e is d e PJ P \p is ~P) 

Furthermore 

1 /° 1 " \ 

A=- 10 p I C~\ (5.21) 

^ \is p 2 pvj 

where we have used that C is symmetric. In particular, we confirm the general relations 


AC = CA T , «Ja, 0 ; 9a',o)) = (AC)aa' , (5.22) 

compare with Appendix If of [14]. 

For the second order expansion we start from the chain rule, 

(dpf\ fdpP dpis d p p\ (dpf\ 

( d r f I = I d r p d r u d r p \d u f . (5.23) 

\9ef J \deLI O e lS d e P J \dpf J 

The Jacobi matrix on the right can be obtained from Eq. (5.20), and together with Eqs. (5.10) 
and (5.13), one arrives at 


fd p (fo)\ /((/o;Po))\ 

U(/o> ]=C- 1 ((/ 0 ;r 0 )) 
\de(fo)J \((/oi e o))/ 

In index notation, 


fdp(fo',ho)\ / ((fo]h 0 ;po))\ 

and d r {f 0 ]ho) = C' 1 ({fo;ho;r 0 )) 

\de{fo] h 0 ) J \((/o; h 0 ; e 0 »/ 


(5.24) 


3 


<(k) = E* 6 " 1 )™' «/o;««'.»» 

a '=1 


3 

and d Ua (f 0 ; h 0 ) = ^ (C~ 1 ) aa > «/o; ho; 9a', o)> • (5.25) 

a'=1 
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This relation allows us to obtain second derivatives as 

3 

du a d u ^(fo) = ^ ( C '“ 1 )aa'(C , “ 1 ) ry ' (((/o;3a',o;5V,0>)) . (5.26) 

a', y=i 

In particular, according to Eqs. (5.14) - (5.16), we obtain the Hessians of the currents 

= (C- 1 H“ C- 1 )^,. (5.27) 

As required, the matrices on the right are symmetric, since C is symmetric. 

A matrix R is introduced in [14] for the transformation to normal modes. Specifically, 
R diagonalizes A and satisfies 

RAR~ 1 = diagonal, RCR T = t. (5.28) 

Inserting Eq. (5.27) into the definition of the coupling matrices, 

3 

G a = IJ2 RoqJ R^H^RT 1 , (5.29) 

a'=1 

and using that C~ 1 R~ 1 = R T leads to 

3 

g° = Raa ' R na ' RT • ( 5 - 30 ) 

af =1 


Special case v = 0. The next step is to take u — 0, which implies r = (rj) = 0, 
(ry p 0 ) = 0, and (rj ; e 0 ) = 0. The inverse matrix of C simplifies to 


c =r 


(( e oUo)) 0 — ((p 0 ;e o)) 

0 0 0 


0 0 0 N 
0 1 0 


«Po; eo» o ;(/)„: p a )) I ^ r °’ r °^ \0 0 0 


(5.31) 


with 


r — ((po; Po))(( e o! e o)) — ((Po! e 0 ))‘ 


Note that T is invariant under e 0 —> e 0 — ppo- 
Using Eq. (5.21) for A at v = 0 results in 


0 


0 


0 


A = — ((e 0 - pp 0 ; e 0 )> 0 -((p 0 ; e 0 - pp 0 )) | + 

'00 0 


'0 1 0 N 
/3((r 0 ;r 0 )) 1° “ o. 


The eigenvalues of A are (—c, 0, c) with c the adiabatic sound speed, 


(5.32) 


(5.33) 


c = ^( r (( r oUo))) 1 / 2 ((e 0 -pPo;e 0 -ppo )) 1/2 


(5.34) 
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Through the relations (5.28) the matrices A and C determine R up to a global sign. 
One obtains 

'(i’-il 
R = I {i’ol 


(5.35) 


(Vh 

where the notation (• | denotes a row vector, and 

'~P" 


tpo = ((eo - lip 0 ;e 0 -ppo))) 1/2 ( 0 


Vv = ((e 0 -ppo\e Q - ppo)) 1 (2T) 


-1/2 /o r \ —1/2 


((e 0 — ppo] e 0 )) 

0 

,-((po; eo — PPo)), 


+ (2((ro;r 0 )))) 


- 1/2 


'O' 

V (5.36) 

for a = ±1. The corresponding inverse matrix, containing the eigenvectors of A as columns, 
is given by 

R- 1 = (iV’-r) |Vh) |Vh» (5.37) 

with 

((do; eo — PPo)}' 

0 

((e 0 — ppo] eo)) t 

Vv = ((eo - P po] e 0 - p po))~ 1/2 (^r ) i/2 | 0 ] + (| ((r 0 ;r 0 ))) 

\dy 


4>o = ((eo - P Po] e 0 - p po)) 1/2 


T/2 


'O' 


(5.38) 


Finally, the G coupling matrices can be determined by inserting i? into Eq. (5.30). First, 
we calculate the entries of R'H P RA , 

$o|Wo> = 0, 

(-0o|?7 p |Vv) = /T 2 (2 ((e 0 - ppo] e 0 - PPo)) ((r 0 ; r o)))~ 1/2 cr , (5.39) 

(Vv|77 p |Vv) = -/3~ 2 (((e 0 - ppo] e 0 - ppo)) {(ro\r 0 )) T) 1/2 ((p 0 ; e 0 -ppo)) <W & 
for cr, cr' = ±1. The entries of RRJRA are 

(VblWo> = o, 

(^o|H r |4)=r 2 (2r)- 1 / 2 , (5.40) 

(i/v|?G|i/v) = -/r 2 r -1 ((p 0 ;e 0 - ppo)) 


and the entries of Rl-L e R T are 

(Am^o) = o, 

(Vh|77 e |-i/v) = /3 -2 (2 ((e 0 -pp 0 ;e 0 - pp 0 )) ((r 0 ; r 0 )))“ 1/2 per , 

(VblT^lVv) = /3~ 2 (((e 0 - ppo; e 0 -ppo)) ((ro\r 0 )) T) 1/2 ((e 0 -2 pp 0 ; e 0 -ppo)) <W a • 

(5.41) 
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Using Eq. (5.30) and (5.34), one obtains for the coupling matrices 


G° 


(-1 0 0 \ 

^-((e 0 -/ipo;e 0 -^po))“ 1/2 0 0 0 

Zp V 0 0 1/ 


and 


G 1 


c 

w 


(( e 0 — PPol e 0 — PPo)) 


1/2 




0 1 \ /0 0 0 \\ 

0 0 + 0 0 1 , 

0 3/ \0 1 0/ / 


where 

T = ~{(pa\e a - ^Po»(2r) _1/2 . 


(5.42) 


(5.43) 


(5.44) 


The matrix G l is determined by G l = — (G 1 ) 7 ”, where T denotes the transpose relative 
to the anti-diagonal. As it has to be, the G coupling matrices are symmetric. 

Note that the thermodynamic averages and cumulants appearing in this paragraph can 
be obtained as appropriate derivatives of the free energy F(p, u, /3) with respect to p, v and 
/3, see Eqs. (5.4) and (5.9). For example, evaluated at v = 0 


((e 0 - ppoi e 0 - ppo)) = -8% (/3F ), (5.45) 

((p 0 ; e 0 - ppo)) = dpd^F , (5.46) 

((Po; Po)) ((e 0 ; e 0 )) - ((p 0 ; e 0 )) 2 = P~ l d}((3F) tf^F - (dpd^F) 2 . (5.47) 

Also expressions involving r 0 , as for example ((ro',r 0 )), can be written as derivative of F 

with respect to v evaluated at v — 0. The numerical computation of F(fx,u,/3) will be 
described below. 

We have arrived at two important qualitative results. Firstly, Gg 0 = 0 always and, 
secondly G\ v > 0, at least for low temperatures, as to be discussed below, and most likely 
in the regime characterized by Eq. (4.8). 


Asymptotic scaling. We define the low temperature correlator by 

-Sit {j,t) = (|p/(t),G(t),e J (t));(p o (0),r o (0),e 0 (0)|)^ (5.48) 

as a 3 x 3 matrix, denoting by |-) the column vector and by (-| its transpose. The average 
is with respect to the thermal state (5.1) at v = 0 and the dynamics is governed by H\ t , 
compare with (4.7). We want to relate S\ t to the DNLS correlator defined by 

= {\Pj(t),fj(t), ej(t))\ (p 0 (0), r 0 (0), e 0 (0)|) M ^ • (5.49) 

Now the dynamics is generated by the DNLS hamiltonian H, compare with (2.16). In (3.2) 
we only considered the correlator for (p 3 , e 3 ). But now we include also the phase difference 
fj = Q(<Pj- i-i — <Pj). We claim that, under the condition (4.8), in approximation 

S(j,t) (5.50) 
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The asymptotics of the correlator on the right can be obtained from nonlinear fluctuating 
hydrodynamics. We use the matrix R , see (5.35), to transform to normal modes and define 


t) = R t) R t , S\j, t) = R S(j, t) R t . (5.51) 

Sf t (j, t) is approximately diagonal with matrix elements 

•S'lt.aa'C 3, t) ^ 5 aa ,f a (j, t ). (5.52) 

The scaling form of ) is explained in [14, 22], For the sound peak one finds 

f a (x,t) ~ (X s t)~ 2/3 f K pz{(X s t)~ 2/3 (x - act)) , (5.53) 

a = ±1, with the non-universal scaling coefficient 

A s = 2v / 2 \Gf a \ . (5.54) 


The universal scaling function /kpz is tabulated in [8], denoted there by /. /kpz > 0, 
/ dx f K pz(x) = 1, f K pz(x) = / K pz(-x), f dxx 2 f KPZ (x ) ~ 0.510523. /kpz looks like a 
Gaussian with a large |x| decay as exp[—0.295|a;| 3 ]. Plots are provided in [23, 8]. The shape 
function for the heat peak is more easily written in Fourier space, 

/o(M) - e _|fc|5/3Ahi , (5.55) 


with the non-universal coefficient 

poo p 

Ah = A/ 2/3 (G° ct ) 2 (47r) 2 / dtt~ 2/3 cos(2vrct) / dx f KPZ (x) 2 . 

Jo Jr 

Numerical evaluation of the free energy. To compute the partition function 

„ N -1 


(5.56) 


v,P) = / p d%d Pj , 


(5.57) 


7=0 


we first switch to polar coordinates (pj,<Pj)- For the case u = 0 we explicitly calculate the 
(Pj integrals [24], i.e., we integrate with respect to fj on [—7r, 7 t]. This leads to 


r N-l 

Z N (fi,0,/3)= / Yl K(p J+ u pj) dp 

J 7=0 


(5.58) 


with the transfer operator or kernel K(x,y) = Ki(x,y)Ko(y) and 

Kfay) = 2wI 0 (pX.y/xy) , K 0 (y) = e^H^X 

Io is the modified Bessel function of the first kind. Then 


(5.59) 


lirn i log Z N (n, 0, P) = log (A max (/l)), 


AT—>-oo 


(5.60) 
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where A max denotes the largest eigenvalue. We follow the ideas in [25, 26] and use a Nystrom- 
type discretization for the kernel. Given a Gauss quadrature rule as 



f(p) e P * 9 ( p dp w ^Wif(xj) 

2=1 


(5.61) 


with positive weights uy and base points x\ [27], we construct the symmetric matrix 

{Kl (Xi , X V ) y/Wi w v ) ” . /=1 (5.62) 

and calculate its largest eigenvalue, denoted Ai. Then 

log(A max (/l)) ~ P Iy + log Ai. (5.63) 

Numerically, we observe exponential convergence with respect to the number of quadrature 
points. At P = 15 and p = 1 for example, n = 16 points suffice for double precision 
accuracy. For non-zero v we proceed analogously. The angular integral is no longer given 
by a special function and we have to determine it numerically. 

Based on (5.60), one then obtains thermodynamic averages and cumnlants as appropri¬ 
ate derivatives of the canonical free energy 

F(p, v, 0) = -/T 1 lim A log Z N (p, u,P) (5.64) 

l\l _V^-N 1 ’ 


with respect to p, u, and p. For example, we determine p numerically such that (pj) = 1, 
as summarized in the following table for several values of P: 


p 

l 

2 

5 

10 

15 

20 

100 

200 


1.05627 

1.18426 

1.08815 

1.03863 

1.02489 

1.01839 

1.003562 

1.001774 


v- 



Figure 2: The chemical potential p as a function of P such that (pj) = 1, for the DNLS 
with parameters m — 1 and g — 1. 

Fig. 2 visualizes p as a function of p. 
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Low temperature expansion. For the special case v = 0, we expand the entries of C, 
A, R , G in terms of /3 _1 at fixed value of the average density p. For that purpose, we expand 
the hamiltonian around the ground state by regarding r 3 and Zj = Pj — p as small. We switch 
to the moving frame, which corresponds to replacing H by H = H — gp (U — pN) — \gp 2 N 
with the ground state energy subtracted. Taylor expansion up to third order results in 

N -1 

H = H cx + ^2 ( - 16 m.p°- (%+i + z j)( z J+ 1 _ z i ) 2 + ^( r j) + ■ (5.65) 

j=o 

The expansion hamiltonian H ex is given by 

N—l 

H ex = C V *3 = h V ) (P + U Z 3 + Z 3+ 0) + 8ife(^'+l - Z if + V( Z j) ( 5 ‘ 66 ) 

3=0 

with potential 

V{x) = \gx 2 for x > —p, V(x) — oo for x < —p. (5.67) 

The boundary condition enforces the constraint Zj > —p. Note that the expansion does 
not depend on (p. Clearly H ex is stable and we have found a state of minimal energy. The 
term cubic in Zj is unstable, but will be stabilized by higher order terms in the expansion. 
These are expected to be small corrections to the quantum pressure ^H( z j +i — z j) 2 and 
hence will be neglected. 

The conserved fields of H cx are Zj, r 3 , tj. Eventually we want to relate them to the 
original held variables. z 3 transforms to pj since the constant shift by p drops out from the 
cumulants below. Concerning the local energy, 

ej ~ tj + gpZj + \gp 2 , that is, tj ~ e 3 — gp pj + const, (5.68) 

such that in a cumulant tj approximates ej — p Pj. 

We use H ex to evaluate the partition function 


Zn(P>,P) — 


,. N -1 

/ e -/3(ifex-AEf=7 1 ^) TT dzjdrj. 

Jr™ , A =o 


(5.69) 


The Tj integrals can be obtained in closed form, and for the z 3 integrals we proceed as 


d2 ym f(z) “ J s dt ySn W°>+/'(°) * +1/» * 2 + • • ■) = m + § /"(o> +.... 

(5.70) 


One obtains 

( z j) 

Setting (zj) = 0 then yields 


I-id obrnor). 

A = (2po)- i + o<r 2 ). 


(5.71) 

(5.72) 
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To leading order the average energy is given by 


fe> = S + (i - |(f+p) _1 )r 1 + °(^ 2 ) ■ ( 5 - 73 ) 

and inserting (5.72) results in 

(e,) = /3 _1 + 0(/3~ 2 ) for (zj) = 0 . (5.74) 

Using relations analogous to (5.9), one obtains 

(( z o', z o)) = ^ (3 1 + 2^(p+p) 13 2 + 0(/3 3 ), (5.75) 

{{to', *o)) = g/3 1_ ^(g+P) P 2 + 0((3 3 ), (5.76) 

(( c oUo )) = ^~I3 X + §(l + P 2 (f + p) ^j(3 2 + 0(13 3 ), (5.77) 

and inserting (5.72) results in 

(( z o]Zo)) = + l(gp(3)~ 2 + 0(/3 ~ 3 ), (5.78) 

((-o; e 0 )) = O + 0(/T 3 ), (5.79) 

«Co; e 0 )) = (3 2 + 0(/3 3 ) (5.80) 

for (zj) = 0. It follows that 

((e 0 - gz 0 ;c Q -ft zq)) = /3~ 2 + 0(/3~ 3 ) (5.81) 

and 

r = 1 /T 3 + 0(p~ 4 ) • (5.82) 

The variance of the stretch is 


((r 0 ; r o» = ™(§ + p) 1 P 1 + 0(/3 2 ), ((r 0 ; r 0 )) = jf3 1 + C?(/3 2 ) for (^) = 0 . (5.83) 

We have collected all ingredients to calculate the remaining thermodynamic quantities. One 
obtains for the square of the sound speed 

° 2 = ^(fr + gp) + 0(P~ 1 ), c 2 = 3-gp + 0(r 1 ) for (zj) = 0 . (5.84) 

The coupling matrix G° is 


G° 


\^Ui3 + 9p) 



+ 0((3 


-i\ 


(5.85) 


and the inner term of G 1 defined in (5.44) is 

T = K? + p)~\‘2gf3)- 1/2 + 0(r 3/2 ) ■ (5.86) 

In particular, > 0 as asserted above. 
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The square of the matrix elements of R 1 are the Landau-Placzek ratios. For the central 
column of R one obtains 


V>o = (~\{ji + gp) 1 ,0,1 - 1(p, + gp) X ) T ^ 1 + 0(fi 2 ), 
%= ( - jjp, 0, l) T /3 _1 + ( 9 (r 2 ) for (zj) = 0 , 

and for the left and right column 

i,„= ((2 9 /?r ‘/ 2 + 0(/r 3/2 )) (o) +((i+p)" 1 / 2 (2/S/m )- 1 / 2 + 0(r 3/2 ) 




(5.87) 


'O' 

u ] , (5.88) 


( (2gf3)~ 1 / 2 + C>(^-3/ 2 ) \ 


Vv = 


for (^) = 0 . 


(5.89) 


(2pP/m)- l / 2 o + 0(p-V 2 ) 

\i(s 9 )- i/2 r 3/2 + o(r 5/2 )/ 

In particular our results imply that, to leading order in /3 _1 , the peak weights for the p-p 
correlations are given by 

((2g/3)- 1 ,(2gp/3)- 2 ,(2gl3)- 1 ) , (5.90) 

for the r-r correlations by 

(2p / d/m)- 1 (l,0,l), (5.91) 

and for the e-e correlations by 

((8 gp 2 )- 1 ?- 3 , r 2 , (8 gp 2 )- 1 ^) • (5-92) 


In the simulations [6] the p-p correlations showed no heat peak, which came unexpected 
at first glance. From our low temperature expansion one concludes that the heat peak is 
suppressed by a relative factor of 1//3, 1 //3 — 0.005 in [6]. 


6 Molecular dynamics simulations 

Time evolution. To solve the DNLS differential equation (2.3) numerically, the symplec- 
tic fourth-order symmetric Runge-Kutta-Nystrom method SRKNg by Blanes and Moan [28] 
is found to be very useful, see also [29]. It can be regarded as generalization of the Strang 
splitting technique. Specifically, a canonical approach for solving the DNLS is to split the 
hamiltonian (2.1) into a kinetic and nonlinear part, H = T + U, with 

N -1 N -1 

i 2 ’ (6i) 

3=0 3=0 

Then the flow over time t (separately induced by T and U) is exactly given by 

: 4> k ^ e - i(1 -™s{2nk/N))t/m^ k ; ( 6 . 2 ) 

(6.3) 
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where denotes the discrete Fourier coefficient k of ('0i)j=o,...,JV-i- The symmetric 

Runge-Kutta-Nystrom method approximates a time step h by the composition 

° Kh ° ' ■ ■ ° ° W ° <„ (6.4) 

with appropriate coefficients a*, 6j. They satisfy the symmetry conditions = a* and 

fcs + 2 -i = b t . In our case the number of stages s = 6, and the numerical values of a* and 
bi for SRKNg are tabulated in [28]. In the simulations below the system size (number of 
lattice sites) is always N = 4096, and the time step h=\. 


Thermodynamic equilibration. For each simulation run, we draw an initial state from 
the canonical ensemble 

JV-l 

Z N (n, P)~ l d pj dqj (6.5) 

3=0 

for a given inverse temperature f3 and chemical potential /i, where H is the DNLS hamil- 
tonian (2.10). To obtain such a state, we discretize the following fictitious overdamped 
Langevin dynamics [30] (also known as biased random-walk) 


dpj(r) = ~dpj(H - /iN) dr + d W P}j (r ), 

d qj(r) = -d qj (H - /iN) dr + ^ dl V q ,j{r) , j = 0, ..., N - 1 


( 6 . 6 ) 


where W p j(t ) t >o and W q j(r) T > o are standard Wiener processes. Numerically, we use the 
Euler-Maruyama method with step size At, 


pT 1 = p]- - "M) + • 

9J+ 1 = i] - Ar d tl (H - M N) + GJj , 

where the • and G v qj are independent standard Gaussian variables. In our implementa¬ 
tion, we set At = T and perform 1024 such steps. We start with a sample drawn from the 
canonical ensemble of the quadratic hamiltonian H- 2 , which is obtained from H ex of (5.66) 
by dropping the cubic term r?(zj + Zj + 1 ). H 2 is diagonalized by Fourier transformation. 


Low temperature dynamics. Numerically we simulate the original DNLS (2.3), with 
initial states drawn from the canonical ensemble of H, and determine S\j,t), compare 
with (5.49) and (5.51). We average over 10 6 simulation runs. Due to the periodic boundary 
conditions, the sum of the stretches fj is necessarily zero or a multiple of 27T. At low 

temperatures there are no umklapp processes, the sum is zero, meaning that with respect 
to fj one simulates a microcanonical ensemble. Since in our theory the canonical average 
is used, we shift the numerical Sff(j,t) a posteriori by the small value ^((ro;ro)). 

Following the proposal in [19], the violation of the conservation law for fj could be 
measured quantitatively by considering 

N- 1 JV-1 

T(X) = ^f(j^) = 53<r,-(t);f°( 0 )). (6.8) 

3=0 j =0 
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At low temperatures, f 3 = r 3 and T(f) = 0, for practical purposes. However at high 
temperatures T(f) increases and reaches its canonical average value for large t. A similar 
consistency check is based on the matrix sum rule 

N -1 

= (6.9) 

j =o 

which holds only if the fields are conserved. For the numerical simulations with j3 — 15 
(see below), the difference is ~ 0.03 with respect to the 3x3 matrix 2-norm at all recorded 
times t. 


•sfiU.t) 



t=0 


t=256 

- t=512 

- t=1024 

- t= 1536 


Figure 3: Equilibrium two-point correlations <5,, (j, t) of the discrete NLS with m = 1 and 
g — 1, showing the right-moving sound peak at different time points. Initial states have 
been drawn from the canonical ensemble (3.1) with parameters f3 = 15 and g = 1.02489. 
The gray vertical lines show the theoretically predicted sound speed c, and the dashed red 
lines are rescaled KPZ functions. 

Fig. 3 visualizes the right-moving sound peak of Aj, (j, t) at different time points, with 
initial states drawn from the canonical ensemble of H with parameters (3 = 15 and g = 
1.02489. Following (5.53), the red dashed lines are KPZ functions scaled as 

(A s t) _2/3 /KPz((A s t) _2/3 (x - ct )) (6.10) 

with c = 0.9556 the theoretical sound speed determined by (5.34). The nonuniversal 
coefficients A s have been fitted to minimize the expression 

N -1 

l^nC ~ (Ast)“ 2/3 /KPz((A s t) _2/3 (j - ct)) | . (6.11) 

l=o 

At the optimal value of A s the error in (6.11) is approximately 0.1 and the largest term 
in the sum (6.11) is approximately 2.5 x 10 -4 at t = 1536, thus confirming the scaling 
exponent and the shape function. The optimal values are recorded in the following table: 
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t 

256 

512 

1024 

1536 

theory 

A s 

0.921 

1.083 

1.186 

1.190 

0.2846 


Up to the largest accessible time, the value seems to stabilize. However, the theoretical 
KPZ scaling prediction 

\ s = 2V2\G° a \, (6.12) 

shown as last entry in the table, still deviates by a factor 4. Such discrepancy, almost 
perfect scaling plot and substantial deviation for the non-universal coefficients, has been 
noted before. We refer to [31] for a discussion. 



Figure 4: Heat mode Sl 0 (j,t) of the discrete NLS, for the same numerical simulations as 
in Fig. 3 at /3 = 15. The dashed orange curves show the fitted Levy distribution in (6.13). 


The (relatively broad) heat mode for the same simulations is visualized in Fig. 4. One 
notices some feedback from the sound modes, which seems to weaken as time progresses. 
According to (5.55), mode coupling predicts 

Si 0 (j,t) = (A h t)- 3/5 /L )5 /3((A h i)- 3/5 j) (6.13) 

with /l iQ the symmetric a-stable distribution, also known as ct-Levy distribution, and shown 
as dashed orange lines in Fig. 4. The coefficients Ah are obtained by minimizing 


N -1 


E kooOU) - (Ah«)- 3/5 /L,5/ 3 ((Ah«)- 3/5 j) 

3=0 


(6.14) 


with the result 
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To probe the domain of validity for nonlinear fluctuating hydrodynamics we move to 
even lower temperatures and set /3 = 200, still at ( pj) = 1, which are the same parameters 
as in [6]. In the top row of Fig. 5 we display the density-density correlation. The central 
peak has dissolved into a bump with rapid oscillations and weakly subballistic spreading. 
The actual heat peak would have a relative amplitude smaller by a factor of 0.005 and is 
thus not visible. The two sound peaks are well separated from the central bump. In [6] their 
structure has been studied in detail for sizes roughly double compared to the one used for 
Fig. 5. There, k , u j space is used, which obstructs a direct comparison. But the numerical 
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Figure 5: Equilibrium density correlations S pp (j,t ) (top row), sound mode (ji, t) (middle 
row) and heat mode SnJj, t) (bottom row) of the DNLS at canonical ensemble parameters 
P = 200 and /i = 1.001774. 


evidence [7] points towards quantitative validity of KPZ scaling. Apparently each sound 
peak is still governed by a scalar noisy Burgers equation. The energy-energy correlations 
look very similar to the one in the top row of Fig. 5. 

Also other linear combinations of the conserved fields can be considered. In the middle 
row of Fig. 5 we display S\ x . As expected only the right sound peak is left, but the central 
bump remains as before. The bottom row of Fig. 5 visualizes Sqq. Now the central bump 
with rapid oscillations has disappeared, giving space for a structure which might be the 
remnant of the heat peak. On the other side, peaks like the sound peak pop up. Apparently, 
we see structures which are not accounted for by the hydrodynamic fluctuation theory. 

At such ultra-low temperatures, the dynamics is no longer sufficiently chaotic, as im¬ 
plicitly assumed when deriving fluctuating hydrodynamics. For a more precise dynamical 
underpinning the obvious candidate is the harmonic approximation with hamiltonian H 2 , 
compare with last paragraph of Sect. 5. In MD simulations of the dynamics for Ho at the 
same parameters, one also observes structures of a similar form as the central bump. A 
direct comparison seems to be difficult. The best studied example of chain dynamics at 
ultra-low temperatures is the FPU /3-chain, size order 10 3 and times up to 10 6 , however at 
a temperature still a factor 1CU 1 lower than studied here [32], For the FPU /3-chain the 
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harmonic approximation is only a small part of the full story, which cautions us to move to 
rapid conjectures. As a general rule, for integrable systems one is left with a case by case 
study. To provide an additional piece of evidence we study an integrable version of DNLS 
in the following section. 


7 The integrable Ablowitz-Ladik model 

The Ablowitz-Ladik (AL) model [33, 18] is an integrable version of the lattice NLS: 

1 37% = “573 A % + 5 » l%f (%+1 + %-1) ■ 


(7.1) 


ft is well known [18, section 3.4] that (7.1) has a hamiltonian structure with non-canonical 
Poisson brackets. Specifically, the hamiltonian in our variables is 

H>11 = J2 e t with e f = u i + v i+' v i) - M 1 - l 9 'm(u 2 + vf)) (7.2) 

jez 


where we split the wavefunction into its real and imaginary part as 

% = 775 ("1 + %)' 


(7.3) 


Note that the logarithmic term in (7.2) requires \ gm(u 2 + v 2 ) = gm\^j\ 2 < 1. The 
system (7.1) can then be written as 


^ Uj = Cj d v ^\ £ Vj = - Cj d Ui H 




d t ■? 


-j ^ Uj J 


(7.4) 


with 


Cj — 1 — | gm(u 2 + v 2 ). (7.5) 

Due to integrability, the Ablowitz-Ladik model has an infinite number of conservation 
laws. We just mention the density 

pf = l (Vh+i + V’l+i fpj) = K u j +1 uj + v j+l Vj) 

with its corresponding current 

J pg = H£~ h\^\ 2 ) fe- 1 rf+1 - ^j-1 V’j+i) 1 

for which the following conservation law holds, 


d al . ^-al _ ^-al _ ri 

d tPj ^P,i+1 Jn.i u ' 


P,3 


The local quantity 


,al 


WJ = log(cj) = log (l - grn |^j 2 ) 
is also conserved, with current 

= 19 W-i Wj -1 - Vh-i ^;_i) • 

Energy is locally conserved with energy current 

/yal __ 1 _ / /ral , 1 /T'al \ 

— mWp,j ‘ 


(7.6) 

(7.7) 

(7.8) 

(7.9) 

(7.10) 

(7.11) 
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Canonical variables and numerical procedure. We apply the ideas in [29, 34] to 
derive a hamiltonian structure with canonical Poisson brackets. First, one can absorb the 
central oscillation by defining, compare with [34], 

Xi (t) = e^’X*), (7.12) 


which transforms (7.1) into 

1 f t Xj = ( 2m + I 9 \Xj?) (Xj +1 + Xj- 1) • 


We define 


and set 


er (x = 


log(l — x) \ 1/2 


x 


for x G (—oo, 1) 


Qj = y2Re(xj)v{gm\xj\ 2 ), 
Pj = V21m(xj)cr(gm\xj\ 2 )- 

Then (7.13) is equivalent to the hamiltonian system 


(7.13) 

(7.14) 


(7.15) 


(7.16) 


with the hamiltonian 

# al = T (l 9 m (pj + i + q'j+i)) r(l gm(p 2 j + q])) (p j+1 Pj + q j+1 qj ) (7.17) 

jez 

and 

/ Q~ X \ 1 /^ 

t(x) — f-——j = 1 — | x + 0(x 2 ) for i6R. (7-18) 

Note that H a[ is symmetric under pj yy qj. The inverse transform to (7.15) is 

Xj = [qj r(\ gm(p 2 + qj)) + i Pj r(| gm{p 2 + qj]) ]. (7.19) 

We apply the Stormer-Verlet scheme [29] to numerically integrate (7.16), and use New¬ 
ton iterations to solve the resulting nonlinear systems of equations. For evaluating corre¬ 
lation functions we transform back to the physical coordinates 'ipj. Since the procedure is 
computationally more demanding compared to the nonintegrable DNLS, we use a smaller 
system size N = 1024 for the Ablowitz-Ladik model, and correspondingly a smaller maxi¬ 
mum correlation time t = 512. We average over 10 5 simulation runs. 

Different from the DNLS, we equilibrate the system by simulating long Hamiltonian 
trajectories and calculating correlation functions for consecutive, disjoint time intervals of 
the trajectory. We set m = g = 1, and choose an initial state such that the average density 
(pf) ~ \ in order to satisfy the constraint gm^j 1 < 1. 
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Figure 6: Equilibrium time-correlations of the Ablowitz-Ladik model at inverse temperature 
/3 = 15. Three time points are superimposed. The faint gray curves are fitted Gaussian 
functions. 
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Figure 7: Equilibrium time-correlations of the Ablowitz-Ladik model at inverse temperature 
f3 = 200 and two time points t = 256 and t = 512 superimposed. 


MD simulation results. Fig. 6 shows the equilibrium time-correlations of the Ablowitz- 
Ladik model at inverse temperature (3 = 15. The correlations have been shifted by a con¬ 
stant to correct for deviations due to the finite number of samples. As for the nonintegrable 
DNLS, left and right sound peaks are clearly visible. However, the energy correlations have 
a sharp cut-off, and the energy peak structure differs from the one of DNLS. 

Due to iutegrability, one expects ballistic spreading of the correlation functions. As 
quantitative test, we fit Gaussian functions to the right density peaks by minimizing 


N/ 2—1 

\ S pp(3,t) - n P (Apf) _1 /G((A p f) _1 (j - c p t)) | . (7.20) 

3=0 

Here n p is a normalization constant, and the three parameters A p , c p , n p are optimized 
numerically, as shown in the following table: 
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128 

256 

512 


0.0642 

0.0639 

0.0636 

c p 

0.422 

0.421 

0.420 
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The values hardly drift with time, which supports the scaling exponent —1 in (7.20). 

Fig. 7 shows the very low temperature (f3 = 200) equilibrium time-correlations of the 
Ablowitz-Ladik model, with two time points t = 256 and t = 512 superimposed. As for the 
nonintegrable DNLS at (3 — 200, the correlations show quickly oscillating features. 

8 Summary and Conclusions 

Our arguments required some length. So let us first summarize what has been accomplished. 
The DNLS has three dynamical regimes. At density 1 we explore a single temperature in 
each regime. At /3 = 1, and certainly for smaller (3, the dynamics is characterized by 
diffusive transport of density and energy. This is the high temperature regime, which is 
very well confirmed by our numerical simulations. Proceeding towards lower temperatures 
we subdivide into the regimes of low and ultra-low temperatures. 

In the low temperature regime the dynamics is still sufficiently chaotic, but the field of 
phase differences is conserved, although only with high accuracy. The three conservation 
laws have equilibrium time correlations whose structure is identical to the one known from 
generic anharmonic chains. On the theoretical level, the main effort is to obtain the coupling 
coefficients needed as an input to nonlinear fluctuating hydrodynamics. Without these 
coefficients, predictions would be mere guess work. Once the couplings are known one can 
take over the results from [14]. The second order expansion of the DNLS Euler currents is 
more difficult than for anharmonic chains. The reason why we still arrive at fairly explicit 
results, valid over the entire low temperature regime, relies on the very special structure 
of the Euler currents, to recall (J p )^ v ,p = u, {J r )y,,v,p — Ab (Je)n,v,p = /xza The Euler 
currents of anharmonic chains have a very different form. For example J r would be the 
momentum which is itself conserved. Still, from the linearized DNLS Euler currents one 
derives two symmetric sound peaks moving with velocity ±c and a heat peak standing 
still. The heat peak broadens as f 3 / 5 with a Levy | shape function and the sound peaks 
broaden as t 2 / 3 with KPZ shape function. In the low temperature regime, no matter what 
the small initial perturbation, and independent of the temperature, one will always obtain 
a linear combination of these three universal peaks in the long time limit. In our molecular 
dynamics simulations at /3 — 15, the scaling of the sound peaks agrees well with the theory. 
The heat peak data are still noisy, but show already the characteristic slow spatial decay. 
For the density-density correlation function the contribution from the heat peak is a factor 
/3 _1 smaller when compared to the sound peak amplitude. 

In the ultra-low temperature regime the integrable structure of the dynamics becomes 
visible. At [3 — 200 one observes structures which indicate already substantial contributions 
from regular motion. Peaks tend to broaden ballistically. The peak structure depends on 
the particular observable and and on temperature. Universal behavior is lost. As illustrated 
in the top row of Fig. 5, KPZ scaling and integrable correlation structure may coexist. Thus 
to identify the precise border line of the ultra-low temperature regime seems to be more 
difficult than distinguishing between the low and high temperature regimes. 

The next, but challenging, goal is to numerically obtain comparable low temperature 
features for the quantized version of the discrete nonlinear Schrodinger equation, which is 
the Bose-Hubbard hamiltonian. 
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A Speed of sound, R matrix and G couplings 

We record the theoretical c, R , and G for the DNLS at our simulation parameters, where 
we set m — 1, g — 1. The matrix G is specified by the relation G _1 = — (G 1 ) 7 ”, with T 
denoting the transpose relative to the anti-diagonal. The entries are rounded to four digits 
for visual clarity. 

At /3 = 15 and p = 1.02489, the speed of sound is c = 0.9556 and 
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as well as 
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At /3 = 200 and // = 1.001774, the speed of sound is c = 0.9970 and 
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B Average currents 

We establish the relations (5.5). To repeat, the average currents for H\ t differ from the ones 
for H, the latter having zero average always. Since our expressions for the currents contain 
derivatives, we use the smooth version of U, V. Then the Boltzmann factor exp[—/3(iTi t — 


1 Pi ~ r j)\ vanishes at pj = 0 and r,j = ±7r. However the Boltzmann factor 

no longer agrees with the one where H\ t is replaced by H. 

For the density current we obtain 


yN -1 


{ a/ Pj-lPj U ( r j~l )) 


= -r 1 ^ 


(^.^e l3Hlt ) (e^ 


0 =*> 


(B.l) 
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where we used partial integration in the last step. Correspondingly for the stretch current, 


r,j ) ^2 \/dj+i/ Pj U(rj) + 2 \JPj~\IPj U(rj_ i) + V ( Pj 


(<V 


-^it) (e^ p i +f}u ^ r i ) = fj, . 


(B.2) 


For the energy current there are more terms to be considered, 

(J.j) = [/'(>-j_i) + {^( r j-i)u'( r j) + r/'(r,_!)[/(>-,))} 

= -r’V J r'( ft )(a rj _,e-^») + 


— u Z 


— pu. 


-1 


+ W Pj+iM 77 (o) +1 v pj-i/pj u (c-i) 




(B.3) 


Taking the limit to an infinitely high potential step, we conclude that the //^-currents 
satisfy (5.5). The zero //-currents would be obtained through a different limit, in which 
the potential barriers are removed. 
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